using System;
using MathLib;

namespace UraniaLib
{
	/// <summary>
	/// Summary description for UraniaSun.
	/// </summary>
	public class UraniaSun
	{
		public UraniaSun()
		{
			//
			// TODO: Add constructor logic here
			//
		}

		public static double GetSolarMEL(DateTime dDate, bool b0Epoch)
		{
			//Returns the Sun's Mean Ecliptic Longitude for the specified date
			double fLong;
			double fJD;
			double fT;

			fJD = UraniaTime.GetJulianDay(dDate, 0);
			if (b0Epoch == true)
			{
				fJD = fJD - 1.0;
			}
        
			fT = (fJD - 2415020.0) / 36525.0;
			fLong = 279.6966778 + (36000.76892 * fT) + (0.0003025 * fT * fT);
			fLong = Trig.PutIn360Deg(fLong);
			return fLong;
		}

		public static double GetSolarPerigeeLong(DateTime dDate, bool b0Epoch)
		{
			//Returns the Sun's Perigee Longitude for the specified date
			double fLong;
			double fJD;
			double fT;

			fJD = UraniaTime.GetJulianDay(dDate, 0);

			if (b0Epoch == true)
			{
				fJD = fJD - 1.0;
			}

			fT = (fJD - 2415020.0) / 36525.0;
			fLong = 281.2208444 + (1.719175 * fT) + (0.000452778 * fT * fT);
			fLong = Trig.PutIn360Deg(fLong);
			return fLong;
		}

		public static double GetSunEarthEcc(DateTime dDate, bool b0Epoch)
		{
			//Returns the eccentricity of Earth's orbit around the sun for the specified date
			double fEcc;
			double fJD;
			double fT;

			fJD = UraniaTime.GetJulianDay(dDate, 0);
			if (b0Epoch == true)
			{
				fJD = fJD - 1.0;
			}
			fT = (fJD - 2415020.0) / 36525.0;
			fEcc = 0.01675104 - (0.0000418 * fT) - (0.000000126 * fT * fT);
			return fEcc;
		}

		public static double GetEarthObliquity(DateTime dDate, bool b0Epoch)
		{
			//Returns theobliquity of Earth's orbit around the sun for the specified date
			double fObl;
			double fDeltaObl;
			double fJD;
			double fT;

			fJD = UraniaTime.GetJulianDay(dDate, 0);
			if (b0Epoch == true)
			{
				fJD = fJD - 1.0;
			}
        
			fT = (fJD - 2415020.0) / 36525.0;
			fDeltaObl = (46.815 * fT) + (0.0006 * fT * fT) + (0.00181 * fT * fT * fT);
			fObl = 23.43929167 - (fDeltaObl / 3600.0);
			return fObl;
		}

		public static double CalcEccentricAnomaly(double fEGuess, double fMA, double fEcc, double fAcc)
		{
			//Calc Ecctrentric Anomaly to specified accuracy
			double fDelta;
			double fDeltaE;
			double fE;
			double fETmp;
			double fEG;

			fEG = fEGuess;

			fDelta = fEG - (fEcc * Math.Sin(fEG)) - fMA;
			if (Math.Abs(fDelta) > fAcc)
			{
				fDeltaE = (fDelta / (1.0 - (fEcc * Math.Cos(fEG))));
				fETmp = fEG - fDeltaE;
				fE = CalcEccentricAnomaly(fETmp, fMA, fEcc, fAcc);
			}
			else
			{
				fE = fEGuess;
			}
			return fE;
		}

		public static void CalcSunPos(DateTime dDate, DateTime dEpoch, ref double fRA, ref double fDecl)
		{
			double fD;
			double fN;
			double fM;
			double fE;
			//double fTanV2;
			//double fV;
			double fLambda;
			double fBeta;
			double fSolarMEL;
			double fSolarPL;
			double fSunEarthEcc;
			double fOblique;
//			double fAcc;

//			fAcc = 0.0000001;
			fBeta = 0;
			fD = UraniaTime.GetDaysBetween(dDate, dEpoch);
			fSolarMEL = GetSolarMEL(dEpoch, true);
			fSolarPL = GetSolarPerigeeLong(dEpoch, true);
			fSunEarthEcc = GetSunEarthEcc(dEpoch, true);
			fOblique = GetEarthObliquity(dEpoch, true);
			fN = (360.0 / 365.242191) * fD;
			fN = Trig.PutIn360Deg(fN);
			fM = fN + fSolarMEL - fSolarPL;
			fM = Trig.PutIn360Deg(fM);
			fM = Trig.DegToRad(fM);
			fE = (360.0 / Math.PI) * fSunEarthEcc * Math.Sin(fM);
			fLambda = fN + fSolarMEL + fE;
			fLambda = Trig.PutIn360Deg(fLambda);

//			fE = CalcEccentricAnomaly(fM, fM, fSunEarthEcc, fAcc);
//			fTanV2 = Math.Sqrt((1.0 + fSunEarthEcc) / (1.0 - fSunEarthEcc)) * Math.Tan(fE / 2.0);
//			fV = Math.Tan(fTanV2) * 2.0;
//			fV = Trig.RadToDeg(fV);
//			fLambda = fV + fSolarPL;
			UraniaCoord.ConvEclToEqu(fOblique, fLambda, fBeta, ref fRA, ref fDecl);
		}

		public static double CalcSunDistance(DateTime dDate, DateTime dEpoch)
		{
			double fD;
			double fN;
			double fM;
			double fE;
			double fTanV2;
			double fV;
			double fDistance;
			double fSolarMEL;
			double fSolarPL;
			double fSunEarthEcc;
			double fAcc;
			double fOblique;
			double fSMA;

			fAcc = 0.0000001;
			fD = UraniaTime.GetDaysBetween(dDate, dEpoch);
			fSolarMEL = GetSolarMEL(dEpoch, true);
			fSolarPL = GetSolarPerigeeLong(dEpoch, true);
			fSunEarthEcc = GetSunEarthEcc(dEpoch, true);
			fOblique = GetEarthObliquity(dEpoch, true);
			fSMA = 149598500.0;

			fN = (360.0 / 365.242191) * fD;
			fN = Trig.PutIn360Deg(fN);
			fM = fN + fSolarMEL - fSolarPL;
			fM = Trig.PutIn360Deg(fM);
			fM = Trig.DegToRad(fM);
			fE = CalcEccentricAnomaly(fM, fM, fSunEarthEcc, fAcc);
			fTanV2 = Math.Sqrt((1.0 + fSunEarthEcc) / (1.0 - fSunEarthEcc)) * Math.Tan(fE / 2.0);
			fV = Math.Atan(fTanV2) * 2.0;
			fDistance = (fSMA * (1.0 - (fSunEarthEcc * fSunEarthEcc))) / (1.0 + (fSunEarthEcc * Math.Cos(fV)));

			return fDistance;
		}
   
		public static double CalcSunAngDiam(DateTime dDate, DateTime dEpoch)
		{
			double fD;
			double fN;
			double fM;
			double fE;
			double fTanV2;
			double fV;
			double fAngDiam;
			double fSolarMEL;
			double fSolarPL;
			double fSunEarthEcc;
			double fAcc;
			double fOblique;
			double fAngDiam0;
			

			fAngDiam0 = 0.533128;
			fAcc = 0.0000001;
			fD = UraniaTime.GetDaysBetween(dDate, dEpoch);
			fSolarMEL = GetSolarMEL(dEpoch, true);
			fSolarPL = GetSolarPerigeeLong(dEpoch, true);
			fSunEarthEcc = GetSunEarthEcc(dEpoch, true);
			fOblique = GetEarthObliquity(dEpoch, true);
	//		fSMA = 149598500.0;

			fN = (360.0 / 365.242191) * fD;
			fN = Trig.PutIn360Deg(fN);
			fM = fN + fSolarMEL - fSolarPL;
			fM = Trig.PutIn360Deg(fM);
			fM = Trig.DegToRad(fM);
			fE = CalcEccentricAnomaly(fM, fM, fSunEarthEcc, fAcc);
			fTanV2 = Math.Sqrt((1.0 + fSunEarthEcc) / (1.0 - fSunEarthEcc)) * Math.Tan(fE / 2.0);
			fV = Math.Atan(fTanV2) * 2.0;
			fAngDiam = fAngDiam0 * ((1.0 + (fSunEarthEcc * Math.Cos(fV))) / (1.0 - (fSunEarthEcc * fSunEarthEcc)));
			return fAngDiam;
		}

		public static double CalcEOT(DateTime dDate, DateTime dEpoch)
		{
			double fRA = 0;
			double fDecl = 0;
			DateTime dGST = new DateTime(2000, 1, 1, 0, 0, 0);
			double fEOT;
			double fUT;
			DateTime dTmp = new DateTime();
			DateTime dUT = new DateTime();

			CalcSunPos(dDate, dEpoch, ref fRA, ref fDecl);

			dGST = UraniaTime.ConvDecTUraniaTime(fRA);
			dUT = UraniaTime.CalcUTFromGSTMod(dGST, dDate);
			dTmp = dTmp.Add(dUT.TimeOfDay);
			fUT = UraniaTime.ConvTimeToDec(dTmp);
			fUT = 12 - fUT;
			fEOT = fUT;
			return fEOT;
		}

		public static double CalcSolarElong(DateTime dDate, DateTime dEpoch, double fPRA, double fPDecl)
		{
			double fSRA = 0;
			double fSDecl = 0;
			double fE;

			CalcSunPos(dDate, dEpoch, ref fSRA, ref fSDecl);

			fPRA = Trig.DegToRad(fPRA * 15.0);
			fPDecl = Trig.DegToRad(fPDecl);
			fSRA = Trig.DegToRad(fSRA * 15.0);
			fSDecl = Trig.DegToRad(fSDecl);
			fE = Math.Acos((Math.Sin(fPDecl) * Math.Sin(fSDecl)) + (Math.Cos(fPRA - fSRA) * Math.Cos(fPDecl) * Math.Cos(fSDecl)));
			fE = Trig.RadToDeg(fE);
			return fE;
		}

		public static double CalcCRN(DateTime dDate)
		{
			double fCRN;
			double fJulian;

			fJulian = UraniaTime.GetJulianDay(dDate, 0);
			fCRN = 1690 + ((fJulian - 2444235.34)/27.2753);
			return fCRN;
		}
	
	}
}
